# LIBRARIES
rm(list= ls())
gc()

ipak <- function(pkg) {
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

Packages <- c("rio", "haven", "janitor", "stargazer", "foreign", "magrittr",
              "nnet", "MASS", "MNLpred", "ggplot2", "ggplot2", "extrafont", "scales", "effects", "ggpubr",
              "jtools", "DAMisc", "simEd")

ipak(Packages)


set.seed(260423)

#LOAD DATA
setwd(" ") #Insert working directory
load("~/ITANES2018_ready.RData") #Insert working directory
ITANES2018<-ITANES2018_ready

#MODELS
#Model 1 --> Table 1
m1<-glm(turnout~fear_all_dummy_numeric+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+eco_hardship_dummy, data=ITANES2018, family="binomial")
m1.coef= exp(coef(m1))


#Model 2 --> Table 1
m2<-multinom(voto_challenger~fear_all_dummy_numeric+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+eco_hardship_dummy, data=ITANES2018, Hess=TRUE)
m2.rrr = exp(coef(m2)) 

#Model 3 --> Table 1 
m3<-multinom(voto_mainstream~fear_all_dummy_numeric+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+eco_hardship_dummy, data=ITANES2018, Hess=TRUE)
m3.rrr = exp(coef(m3)) 

#Model 4 --> Table A1 
m4<-multinom(voto~fear_all_dummy_numeric+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+eco_hardship_dummy, data=ITANES2018, Hess=TRUE)
m4.rrr = exp(coef(m4)) 

#Model 5 --> Table S1 
m5<-multinom(voto~employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+eco_hardship_dummy, data=ITANES2018, Hess=TRUE)
m5.rrr = exp(coef(m5)) 

#Model 6 --> Table S1
m6<-multinom(voto~fear_all_dummy_numeric*employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+eco_hardship_dummy, data=ITANES2018, Hess=TRUE)
m6.rrr = exp(coef(m6)) 

#Model 7 --> Table 2
m7<-glm(turnout~fear_all_dummy_numeric*eco_hardship_dummy+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent, data=ITANES2018, family="binomial")
m7.coef= exp(coef(m7))

#Model 8 --> Table 2
m8<-multinom(voto~fear_all_dummy_numeric*eco_hardship_dummy+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent, data=ITANES2018, Hess=TRUE)
m8.rrr = exp(coef(m8)) 

#Model 2b --> Table S2 
m2b<-multinom(voto_Lega_challenger~fear_all_dummy_numeric+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+eco_hardship_dummy, data=ITANES2018, Hess=TRUE)
m2b.rrr = exp(coef(m2b)) 

#Model 3b --> Table S2
m3b<-multinom(voto_Lega_mainstream~fear_all_dummy_numeric+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+eco_hardship_dummy, data=ITANES2018, Hess=TRUE)
m3b.rrr = exp(coef(m3b)) 

#Model 4b --> Table S3
m4b<-multinom(voto_Lega~fear_all_dummy_numeric+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+eco_hardship_dummy, data=ITANES2018, Hess=TRUE)
m4b.rrr = exp(coef(m4b)) 

#Model 2c --> Table S4
m2c<-multinom(voto_challenger~fear_activeworkers_dummy+eco_hardship_dummy+employment_status_nounemp+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent, data=ITANES2018, Hess=TRUE)
m2c.rrr = exp(coef(m2c)) 

#Model 3c --> Table S4
m3c<-multinom(voto_mainstream~fear_activeworkers_dummy+eco_hardship_dummy+employment_status_nounemp+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent, data=ITANES2018, Hess=TRUE)
m3c.rrr = exp(coef(m3c)) 

#Model 4c --> Table S5
m4c<-multinom(voto~fear_activeworkers_dummy+eco_hardship_dummy+employment_status_nounemp+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent, data=ITANES2018, Hess=TRUE)
m4c.rrr = exp(coef(m4c)) 

#Model 1d --> Table S6
m1d<-glm(turnout~fear_all_dummy_numeric+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+voto_abstain_2013+eco_hardship_dummy, data=ITANES2018, family="binomial")
m1d.coef= exp(coef(m1d)) 

#Model 2d  --> Table S6
m2d<-multinom(voto_challenger~fear_all_dummy_numeric+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+voto_abstain_2013+eco_hardship_dummy, data=ITANES2018, Hess=TRUE)
m2d.rrr = exp(coef(m2d)) 

#Model 3d --> Table S6
m3d<-multinom(voto_mainstream~fear_all_dummy_numeric+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+voto_abstain_2013+eco_hardship_dummy, data=ITANES2018, Hess=TRUE)
m3d.rrr = exp(coef(m3d)) 

#Model 4d --> Table S7
m4d<-multinom(voto~fear_all_dummy_numeric+employment_status+age+gender+education+residence+union+ideology+euro+immigration+trust+populism+incumbent+voto_abstain_2013+eco_hardship_dummy, data=ITANES2018, Hess=TRUE)
m4d.rrr = exp(coef(m4d)) 



#PRINT PLOTS
#TABLE 1: m1,m2,m3 
stargazer(m1,m2,m3, coef = list(m1.coef,m2.rrr,m3.rrr),  p.auto=FALSE, type="text")

#TABLE 2: m7, m8 
stargazer(m7,m8, coef = list(m7.coef,m8.rrr),  p.auto=FALSE, type="text")

#TABLE A1: m4 
stargazer(m4, coef = list(m4.rrr),  p.auto=FALSE, type="text")

#TABLE A2: m5 
stargazer(m5, coef = list(m5.rrr),  p.auto=FALSE, type="text")

#TABLE A3: m6 
stargazer(m6, coef = list(m6.rrr),  p.auto=FALSE, type="text")

#TABLE S1: m2b, m3b 
stargazer(m2b,m3b, coef = list(m2b.rrr,m3b.rrr),  p.auto=FALSE, type="text")

#TABLE S2: m4b 
stargazer(m4b, coef = list(m4b.rrr),  p.auto=FALSE, type="text")

#TABLE S3: m2c, m3c
stargazer(m2c,m3c, coef = list(m2c.rrr,m3c.rrr),  p.auto=FALSE, type="text")

#TABLE S4: m4c
stargazer(m4c, coef = list(m4c.rrr),  p.auto=FALSE, type="text")

#TABLE S5: m1d,m2d,m3d 
stargazer(m1d,m2d,m3d, coef = list(m1d.coef,m2d.rrr,m3d.rrr),  p.auto=FALSE, type="text")

#TABLE S6: m4d 
stargazer(m4d, coef = list(m4d.rrr),  p.auto=FALSE, type="text")


#FIGURES 

#Figure A1
multi1<-multinom(voto~fear_all_dummy_numeric+eco_hardship_dummy+permanently_employed+atypically_employed+unemployed+age+gender+education+Center+South+union+left+centre_left+centre+centre_right+right+euro+immigration+trust+populism+incumbent, data=ITANES2018, Hess=TRUE)
set.seed(260423)
pred_perceived <- mnl_pred_ova(model = multi1,
                               data = ITANES2018,
                               x = "fear_all_dummy_numeric",
                               by = 1,
                               seed = 260423,
                               nsim = 2000, # faster
                               probs = c(0.025, 0.975))

pred_perceived$plotdata %>% head()
pred_perceived$plotdata$fear_all_dummy<-pred_perceived$plotdata$fear_all_dummy_numeric

pred_perceived$plotdata$fear_all_dummy<-as.factor(as.numeric(pred_perceived$plotdata$fear_all_dummy))
levels(pred_perceived$plotdata$fear_all_dummy)
levels(pred_perceived$plotdata$fear_all_dummy)[1]<-"Not afraid"
levels(pred_perceived$plotdata$fear_all_dummy)[2]<-"Afraid"

png("pp.png", units="in", width=8, height=5, res=600)

ggplot(data = pred_perceived$plotdata, aes(x = fear_all_dummy, shape=fear_all_dummy, 
                                           y = mean,
                                           ymin = lower, ymax = upper)) + 
  geom_pointrange() + # Confidence intervals
  # geom_line() + # Mean
  facet_wrap(.~ voto,  ncol = 3) +
  scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(0, 0.6), breaks=c(0,0.2,0.4, 0.6)) + # % labels
  
  theme_bw() +
  labs(y = "Predicted probabilities",
       x = "")+ 
  facet_wrap(.~ voto,  nrow = 2)+ 
  theme(axis.text.y = element_text( size=14))+ 
  theme(axis.title = element_text(size = 14))+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+ labs(shape = "Fear of job loss")+ theme(legend.text=element_text(size=12))+ theme(legend.title=element_text(size=12))+
  theme(text = element_text(family = "LM Roman 10"))

dev.off()

#Figure 2
set.seed(260423)
fdif_perceived <- mnl_fd2_ova(model = multi1,
                              data = ITANES2018,
                              x = "fear_all_dummy_numeric",
                              value1 = min(ITANES2018$fear_all_dummy_numeric, na.rm=TRUE),
                              value2 = max(ITANES2018$fear_all_dummy_numeric, na.rm=TRUE),
                              nsim = 2000,
                              seed = 260423)

png("fd.png", units="in", width=8, height=5, res=600)
ggplot(fdif_perceived$plotdata_fd, aes(x = categories, 
                                       y = mean,
                                       ymin = lower, ymax = upper)) +
  geom_pointrange() +
  geom_hline(yintercept = 0) +
  scale_y_continuous(labels = percent_format()) +
  theme_bw() +
  labs(y = "Predicted probabilities",
       x = "Vote choice") + theme(axis.text.x = element_text( size=14))+ 
  theme(axis.text.y = element_text( size=14))+ 
  theme(axis.title = element_text(size = 15)) + theme(
    axis.title.x = element_text(hjust=0.5),
    axis.title.y = element_text(hjust=0.5)
  )   +
  theme(text = element_text(family = "LM Roman 10"))
dev.off()



# Figure A2
df1 <- ITANES2018
df1$atypically_employed <- 0
df1$selfemployed <- 0
df1$unemployed <- 0
set.seed(260423)
pred_permanently_employed <- mnl_pred_ova(model = multi1,
                                          data = df1,
                                          x = "permanently_employed",
                                          by = 1,
                                          seed = 260423,
                                          nsim = 2000, # faster
                                          probs = c(0.025, 0.975)) 
set.seed(260423)
fdif_permanently_employed <- mnl_fd2_ova(model = multi1,
                                         data = df1,
                                         x = "permanently_employed",
                                         value1 = min(df1$permanently_employed,na.rm=TRUE),
                                         value2 = max(df1$permanently_employed,na.rm=TRUE),
                                         nsim = 2000,
                                         seed = 260423)


df2 <- ITANES2018
df2$permanently_employed <- 0
df2$selfemployed <- 0
df2$unemployed <- 0
set.seed(260423)
pred_atypically_employed <- mnl_pred_ova(model = multi1,
                                         data = df2,
                                         x = "atypically_employed",
                                         by = 1,
                                         seed = 260423,
                                         nsim = 2000, # faster
                                         probs = c(0.025, 0.975)) 
set.seed(260423)
fdif_atypically_employed <- mnl_fd2_ova(model = multi1,
                                        data = df2,
                                        x = "atypically_employed",
                                        value1 = min(df2$atypically_employed,na.rm=TRUE),
                                        value2 = max(df2$atypically_employed,na.rm=TRUE),
                                        nsim = 2000,
                                        seed = 260423)

df3 <- ITANES2018
df3$permanently_employed <- 0
df3$selfemployed <- 0
df3$atypically_employed <- 0
set.seed(260423)
pred_unemployed <- mnl_pred_ova(model = multi1,
                                data = df3,
                                x = "unemployed",
                                by = 1,
                                seed = 260423,
                                nsim = 2000, # faster
                                probs = c(0.025, 0.975)) 
set.seed(260423)
fdif_unemployed <- mnl_fd2_ova(model = multi1,
                               data = df3,
                               x = "unemployed",
                               value1 = min(df3$unemployed,na.rm=TRUE),
                               value2 = max(df3$unemployed,na.rm=TRUE),
                               nsim = 2000,
                               seed = 260423)

pred_unemployed$plotdata$category <- ifelse(pred_unemployed$plotdata$unemployed == 1, "Unemployed", "Self-employed")
pred_atypically_employed$plotdata$category <- ifelse(pred_atypically_employed$plotdata$atypically_employed == 1, "Atypically employed", "Self-employed")
pred_permanently_employed$plotdata$category <- ifelse(pred_permanently_employed$plotdata$permanently_employed == 1, "Permanently employed", "Self-employed")
library(dplyr)
pred_atypically_employed$plotdata <- pred_atypically_employed$plotdata %>% 
  filter(category != "Self-employed")
pred_permanently_employed$plotdata <- pred_permanently_employed$plotdata %>% 
  filter(category != "Self-employed")

# bind all pps together

pps <- pred_unemployed$plotdata %>%
  dplyr::bind_rows(pred_atypically_employed$plotdata) %>%
  dplyr::bind_rows(pred_permanently_employed$plotdata) 

table(pps$category)
pps$category<-as.factor(as.character(pps$category))


pps$category<-factor(pps$category, levels=c("Unemployed" , "Atypically employed","Permanently employed", "Self-employed"))


png("pp_emp.png", units="in", width=8, height=5, res=600)

ggplot(data = pps, aes(x = category, shape = category,
                       y = mean,
                       ymin = lower, ymax = upper)) +
  #  scale_shape_manual(values = c(18,15,17,20)) +
  geom_pointrange()+ # Confidence intervals
  # geom_line() + # Mean
  facet_wrap(.~voto,  ncol = 2)+
  scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(0, 0.6), breaks=c(0,0.2,0.4,0.6))+ # % labels
  # scale_x_continuous(breaks = c(0:10), limits = c(-0.3, 1.3)) +
  theme_bw() +
  labs(x = "",
       y = "Predicted probabilities")+ scale_color_brewer(type = 'div', palette = 4)+ theme(axis.text.x = element_text( size=14))+ 
  theme(axis.text.y = element_text( size=14))+ 
  theme(axis.title = element_text(size = 15))+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+ labs(shape = "Employment status")+ theme(legend.text=element_text(size=12))+ theme(legend.title=element_text(size=12))+
  theme(text = element_text(family = "LM Roman 10"))

dev.off()


# Predictions 2013 (not included in the paper)
table(ITANES2018$voto_abstain_2013)
ITANES2018$abstain_2013<-ifelse(ITANES2018$voto_abstain_2013=="Abstain", 1, 0)
ITANES2018$Challengerleft_2013<-ifelse(ITANES2018$voto_abstain_2013=="Challenger-left", 1, 0)
ITANES2018$Challengerright_2013<-ifelse(ITANES2018$voto_abstain_2013=="Challenger-right", 1, 0)
ITANES2018$Mainstreamright_2013<-ifelse(ITANES2018$voto_abstain_2013=="Mainstream right", 1, 0)
ITANES2018$Mainstreamleft_2013<-ifelse(ITANES2018$voto_abstain_2013=="Mainstream left", 1, 0)

#Preds perceived from model 4 --> Fig. S1
multi1d<-multinom(voto~fear_all_dummy_numeric+eco_hardship_dummy+permanently_employed+atypically_employed+unemployed+age+gender+education+Center+South+union+left+centre_left+centre+centre_right+right+euro+immigration+trust+populism+incumbent+Mainstreamright_2013+Mainstreamleft_2013+Challengerright_2013+Challengerleft_2013, data=ITANES2018, Hess=TRUE)
set.seed(260423)
pred_perceived <- mnl_pred_ova(model = multi1d,
                               data = ITANES2018,
                               x = "fear_all_dummy_numeric",
                               by = 1,
                               seed = 260423,
                               nsim = 2000, # faster
                               probs = c(0.025, 0.975))

pred_perceived$plotdata %>% head()
pred_perceived$plotdata$fear_all_dummy<-pred_perceived$plotdata$fear_all_dummy_numeric

pred_perceived$plotdata$fear_all_dummy<-as.factor(as.numeric(pred_perceived$plotdata$fear_all_dummy))
levels(pred_perceived$plotdata$fear_all_dummy)
levels(pred_perceived$plotdata$fear_all_dummy)[1]<-"Not afraid"
levels(pred_perceived$plotdata$fear_all_dummy)[2]<-"Afraid"

png("pp2013.png", units="in", width=8, height=5, res=600)

ggplot(data = pred_perceived$plotdata, aes(x = fear_all_dummy, shape=fear_all_dummy, 
                                           y = mean,
                                           ymin = lower, ymax = upper)) + 
  geom_pointrange() + # Confidence intervals
  # geom_line() + # Mean
  facet_wrap(.~ voto,  ncol = 3) +
  scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(0, 0.6), breaks=c(0,0.2,0.4, 0.6)) + # % labels
  
  theme_bw() +
  labs(y = "Predicted probabilities",
       x = "")+ 
  facet_wrap(.~ voto,  nrow = 2)+ 
  theme(axis.text.y = element_text( size=14))+ 
  theme(axis.title = element_text(size = 14))+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+ labs(shape = "Fear of job loss")+ theme(legend.text=element_text(size=12))+ theme(legend.title=element_text(size=12))+
  theme(text = element_text(family = "LM Roman 10"))

dev.off()

# Figure S1
set.seed(260423)
fdif_perceived <- mnl_fd2_ova(model = multi1d,
                              data = ITANES2018,
                              x = "fear_all_dummy_numeric",
                              value1 = min(ITANES2018$fear_all_dummy_numeric, na.rm=TRUE),
                              value2 = max(ITANES2018$fear_all_dummy_numeric, na.rm=TRUE),
                              nsim = 2000,
                              seed = 260423)

png("fd2013.png", units="in", width=8, height=5, res=600)
ggplot(fdif_perceived$plotdata_fd, aes(x = categories, 
                                       y = mean,
                                       ymin = lower, ymax = upper)) +
  geom_pointrange() +
  geom_hline(yintercept = 0) +
  scale_y_continuous(labels = percent_format()) +
  theme_bw() +
  labs(y = "Predicted probabilities",
       x = "Vote choice") + theme(axis.text.x = element_text( size=14))+ 
  theme(axis.text.y = element_text( size=14))+ 
  theme(axis.title = element_text(size = 15)) + theme(
    axis.title.x = element_text(hjust=0.5),
    axis.title.y = element_text(hjust=0.5)   
  )+
  theme(text = element_text(family = "LM Roman 10"))
dev.off()




#Predictions turnout

newdata <- with(ITANES2018, data.frame(fear_all_dummy_numeric = c(0,1), employment_status="Permanently employed", 
                                       age=mean(age, na.rm=TRUE), gender= 1, education=mean(education, na.rm=TRUE),
                                       residence="North",union= 0, ideology="Left", euro=mean(euro, na.rm=TRUE), 
                                       immigration=mean(immigration, na.rm=TRUE), trust=mean(trust, na.rm=TRUE), 
                                       eco_hardship_dummy=mean(eco_hardship, na.rm=TRUE), populism=mean(populism, na.rm=TRUE), 
                                       incumbent=mean(incumbent, na.rm=TRUE)))
set.seed(260423)
preds <- predict(m1, newdata, type="response", se.fit=TRUE)
predf <- preds$fit # predicted
lower <- preds$fit - (1.96*preds$se.fit) # lower bounds
upper <- preds$fit + (1.96*preds$se.fit) # upper bounds
predf


